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Abstract. A numerical method to solve the direct scattering problem for the 
Zakharov-Shabat system associated to the initial value problem for the non¬ 
linear Schrodinger equation is proposed. The method involves the numerical 
solution of Volterra integral systems with structured kernels and the identifica¬ 
tion of coefficients and parameters appearing in monomial-exponential sums. 
Numerical experiments confirm the effectiveness of the proposed technique. 
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1. Introduction 


The problem we are addressing concerns the numerical computation of the scat¬ 
tering data of the Zakharov-Shabat (ZS) system associated to the initial value 
problem (IVP) for the nonlinear Schrodinger (NLS) equation 


(i.i) 


hit + u xx ± 2\u\ 2 u = 0, x £ M, t > 0 
tt(a;,0) = uq(x), 


where i denotes the imaginary unit, u = u(x,t) is the unknown potential, the 
subscripts x and t designate partial derivatives with respect to position and time, 
uq £ L 1 (R) is the initial potential and the ± sign depends on symmetry properties 
of u. The plus sign regards the focusing case and the minus sign the defocusing 
case. 

The solution of the IVP (11.111 can theoretically be obtained by means of the 
so-called Inverse Scattering Transform (1ST) technique [BE]. The 1ST allows one, 
in fact, to obtain the solution of m by means of the following three steps: 

(i) starting from the initial potential uq, solve the Zakharov-Shabat (ZS) sys¬ 
tem associated to the NLS to obtain the initial scattering data; 

(ii) propagate the initial scattering data in time; 

(iii) solve the associated Marchenko equations whose kernels are obtained from 
the initial scattering data evolved in time, to obtain the solution u(x, t) we 
are looking for. 

An effective numerical method to solve steps (ii) and (iii) has been proposed in 
2] under the hypothesis that the initial scattering data are known. In this paper 
we propose a numerical method to solve the direct scattering problem (i) which 
is also of independent interest in some engineering fields M- To the best of our 
knowledge our method is the first numerical method proposed for the computation 
of all scattering data. 

The paper is organized as follows. In Section [2] we recall the ZS system as¬ 
sociated to the IVP for the NLS equation. Then we recall the definition of the 
initial scattering data, i.e. the transmission coefficient T{ A), the reflection coeffi¬ 
cient from the left L( A) and from the right R( A), the bound states {A^} with their 
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multiplicities {irij} and the norming constants from the left {and from the 
right {{T r ) js }. After that, we introduce the initial Marchenko kernels from the left 
idf(a) and from the right Q r (a), the inverse Fourier transform p(a) of i?(A) and the 
Fourier transform t(a) of L( A), respectively. Then, we show that the spectral sums 
from the left St (a) and from the right S r (a) which depend on the bound states 
with the respective multiplicities and norming constants from the left and from the 
right, can be expressed as a difference between the initial Marchenko kernels and 
the inverse Fourier transform p(a) of i?(A) and the Fourier transform £(a) of L( A), 
respectively. As these differences are monomial-exponential sums, their parameters 
and coefficients can be identified by using the numerical method proposed in [HJ (Bi¬ 
section [3] is devoted to the characterization of auxiliary functions which are basic 
to the computation of the initial Marchenko kernels as well as of p(a) and £(a). 
In Section [3] we characterize the scattering matrix, derive and analyze the Volterra 
integral equations of the second kind that characterize the initial Marchenko ker¬ 
nels and formulate the Fredholm integral equations that characterize the Fourier 
transforms p(a) and £(a). The numerical method we propose to obtain the initial 
scattering data is illustrated in Section [5] while in Section [G] we consider two differ¬ 
ent initial potentials for which the numerical results are given in Section 0 Finally, 
we conclude the paper by an Appendix concerning the study of the support of the 
auxiliary functions introduced in Section [3] 


2. Initial scattering data 


Following the 1ST technique, to determine the initial scattering data, we must 
consider the ZS system associated to the NLS m a, that is the system 

(2.1) iJ^(A,s)-V(a)$(A,s)=A$(A,a:), iel 

OX 

where A G C is a spectral parameter and 


( 2 . 2 ) 





with vo = u o in the focusing case and vq = —uo(x)* in the defocusing case. Here 
and in the sequel the asterisk denotes the complex conjugate. 

The initial scattering data are the entries of the so-called scattering matrix and 
the coefficients and parameters of two spectral sums. Denoting by 


S(A) 


(n a ) m\ 

\R{ A) T(X)J 


the scattering matrix, T(A) represents the (initial) transmission coefficient, while 
L( A) and R( A) stand for the initial reflection coefficients from the left and from the 
right, respectively. This matrix satisfies the following symmetry properties ll3j 

(2.3) S f (A)S(A) = S(A)S f (A) =1 


in the defocusing case and 

(2.4) S t (A)JS(A) =S(A)JS t (A) = J 


in the focusing case where I denotes the identity matrix. Here and in the sequel 
the dagger denotes the matrix conjugate transpose. The numerical validity of these 
properties is used in Section [7] to check the effectiveness of our algorithms. 

If T(A) has no poles in the complex upper half plane C + , there are no spectral 
sums to identify. Otherwise, denoting by Ai, ..., A„ the so-called bound states, 
that is the finitely many poles of T{ A) in C + , and by mi, ..., m n the correspond¬ 
ing multiplicities, we have to identify the parameters {n, rrij,Xj} as well as the 
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coefficients {(r^)j s , (r r ) JS } of the initial spectral sums from the left and from the 
right 

n TOjj- 1 s 

(2-5) S e (a) = J2 eiXja 

j—1 s=0 

n rrij -1 8 

(2.6) Sr (a) = Y / e iX * a E 

j=1 s=0 

In (12.511 and (12.61) the coefficients (r^) JS and (r r ) JS are the so-called norming con¬ 
stants from the left and from the right, respectively, and 0! = 1. 

In the 1ST technique, a crucial role is played by the initial Marchenko kernels 
from the left fle(a) and from the right fl r (a), which are connected to the above 
spectral coefficients and spectral sums as follows: 


(2.7) 

fle(a) = p(a ) + Se(a), for a > 0 

(2.8) 

f l r (a) = £(a) + S r (a), for a < 0 

where 

1 r +0 ° 

p(a ) = — / R(\)e iXa d\ = R~ l {R( A)} 

2tt J.oo 

(2.9) 


is the inverse Fourier transform of the reflection coefficient from the right R( A) and 

I r+°° i 

(2.10) 1(a) = — j L(\)e~ iXa d\ = —T {L( A)} , 

apart from the factor 1/27r, is the Fourier transform of the reflection coefficient 
from the left L( A). 

We note that Sl^(a) and fl r (a), respectively, reduce to: 

(a) Se(a) and S r (a ) if the reflection coefficients vanish (reflectionless case); 

(b) p(a) and £(a) if there are no bound states. 


a > 0 

a < 0. 


3. Auxiliary functions 


In this section we introduce four pairs of auxiliary functions and the Volterra 
integral equations that characterize them. Their solution, as shown in the next 
section (see also Bed]), is fundamental for computing the initial Marchenko kernels 
as well as p(a) and £(a). 

Following [7], let us introduce, for y > x, the two pairs of unknown auxiliary 
functions 


K (x,y) 


(K" p (x,y)\ 
\K dn (x,y)J ’ 


K (x,y) 


( K up (x,y)\ 
\K dn (x,y)J ’ 


and, for y < x, the two pairs of unknown auxiliary functions 


M(x, y) 


fM pp (x,y)\ 
\M dn (x,y)J ’ 


M (x,y) 


[M pp (x,y)\ 
\M dn (x,y)J ' 


For the sake of clarity, let us explain how these functions are connected to the 
Jost matrices associated to the ZS system m- 

As in mm, we represent the Jost matrices as the Fourier transforms of the 
auxiliary functions: 
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(3.1) (*(\,x),<S'(\,x)) = e~' XJX + I (K(x, y), K(x, y)) e~ iX3y dy 


(3.2) (<f>(\, X ),H>(\,x)) = e- iXJx + ( M(x,y),M(x,y))e- iXJ ydy , 

J — OO 

from which inverting the Fourier transforms we get 

(3.3) (K(x,y),K(x,y)) = — J [(^(A,*), ¥(A, x)) - e~ iXJ *]e iXJ y dX, 

(3.4) (M(x,y),M(x,y)) = — J [(*(A, x), *(A, x)) - e~ iXJx ]e iXJ « d\. 

Now, for y > x, the pair (A' up , A' dn ) is the solution of the following system of two 
structured Volterra integral equations mm- 


(3.5) 


K up (x,y) + / uq(z) I< dn (z,z + y - x) dz = 0 


/■ 2 (*+?/) 

K dn {x,y )- / v 0 (z) K up (z,x + y- z) dz = \v 0 (\(x + y)) 


while the pair (A' up , AT dn ) is the solution of the system 


r 5 (*+») 

K up (x,y)+ / uo(z)K dn (z, X + y-z)dz = ~±uo(^(x + y)) 

J X 

/»oo 

K dn (x,y) — / vo(z) K up (z, z + y — x) dz = 0. 

< J X 


(3.6) 


Similarly, for y < x the pair ( M up , M dn ) is the solution of the system of two 
structured Volterra equations: 


M" p {x,y)- 


(3.7) 


>(*+i0 


u 0 {z) M dn (z,x + y — z) dz = hu 0 (h(x + y)) 


M dn (x, y) + / v 0 (z)M np (z,z + y-x)dz = 0 

k J — OO 

and the pair (M up , M dn ) is the solution of the following system 

' pX 

M up (x, y)- u 0 (z) M dn (z, z + y-x)dz = 0 

J —OO 


(3.8) 


M dn (x,y)+ ( v 0 (z)M up (z,x + y- z) dz = -\v 0 (\{x + y)). 

J o 0+2/) 


From the computational point of view, on the bisector y = x, it is important 
to note that each auxiliary function is uniquely determined by the initial solution 
or its partial integral energy. In fact, setting y = x in each of the four Volterra 
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systems, we immediately obtain: 


i ,, - , , i r 


(3.9) 

A' dn (x, x) = -v Q {x), 

K up (x, x) = - 

2 J x 

u 0 (z)v 0 (z) dz, 

(3.10) 

K up (x,x) = -^uo(x), 

K dn (x,x) = 

4 . 

/*oo 

/ u 0 (z)v 0 (z)dz, 

' X 

(3.11) 

M dn (x,x) = ~v 0 (x), 

M up (x, x) = 

4 

/ u 0 (z)v 0 (z)dz, 

' —OO 

(3.12) 

M up (x,x) = \u 0 {x), 

M dn (x,x) = - 

-)/ 

>x 

u 0 (z) v 0 (z) dz. 


J — OO 


Moreover, let us mention that the functions K and K, as well as the functions 
M and M, are related to each other. Indeed, in the focusing case the following 
symmetry properties hold true m 


(3.13) 


1 I\ up (x , y) N 
,K dn (x,y). 


'-K d *{x,y)* 
K up (x, y)* 


(M up (x , y 4 
K M dn (x,y). 


M dn (x, y)* 
,—M up (x, y)*, 


while in the defocusing case the following symmetry relations can be proved 


(3.14) 


K up (x, y) 


K dn (x, y)* 


M up (x,y) 


M dn (x,yY 


\IC dn {x,y)J \K up (x, y)*J \M du (x,y)J \M up (x,y)*J 

Remark 1. Let us note that, in virtue of (I3.13I) - (I3.14D . we only need to compute 
numerically systems m and m or systems (EH) and EH and then compute 
the remaining auxiliary functions by resorting to the above symmetry properties . 

Remark 2. If the potentials uq and vq are even functions, the auxiliary functions 
M can easily be obtained from the K functions as follows: 

M up (x, y) = K up (—x, -y) M dn (x, y) = -7? dn (-x, - y) 

M" p (x, y) = —K up (—x, -y) M dn (x, y) = A' dn (-x, -y). 

Similarly, if the potential Uq and vq are odd functions, we have 

M up (x, y ) = K" p (-x,-y) M dn (x, y) = K d »(-x, -y) 

M up (x, y ) = K pp (-x,-y) M dn (x, y ) = A' dn (-x, -y). 

Consequently, in these cases we only need to compute numerically one system, for 
instance system (13.51) . 


4. Initial Marchenko kernels, scattering matrix and Fourier 

TRANSFORMS OF REFLECTION COEFFICIENTS 

This section consists of two parts. In the first part we recall the Volterra integral 
equations that we solve to obtain the initial Marchenko kernels f lt{a) and fl r (a). 
In the second part we explain how compute the scattering matrix and the Fourier 
transforms of the reflection coefficients R and L. 


4.1. Initial Marchenko kernels. Following jT3} 2.50a and 2.50b] we can say that, 
for y > x > 0, the Marchenko kernel f^ is connected to the auxiliary functions A' dn 
and A' dn as follows: 


A' dn (x, z) fle(z + y) dz 


—K An (x, y). 


(4.1) 


&i{x + y) + 
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Similarly, for y < x < 0, the Marchenko kernel fl r is connected to the auxiliary 
functions M dn and M dn in this way: 

(4.2) il r (x + y)+ f M up (x,z)Cl r (z + y)dz =—M up (x,y). 

J — OO 

As a result, assuming known the auxiliary functions, (14.11) and (14.21) can be in¬ 
terpreted as structured Volterra integral equations having the initial Marchenko 
kernels fie and fl r as their unknowns. 

It is important to note that, from the computational point of view, each Marchenko 
kernel can be treated as a function of only one variable, as we only have to deal 
with the sum of the two variables. 


4.2. The scattering matrix and the Fourier transforms of the reflection 
coefficients. Let us begin by recalling that, as proposed in [T3], the coefficients of 
the scattering matrix ^(A) can be represented as follows: 


(4.3) 

(4.4) 

(4.5) 


T(\) = —77T =-7v\ ’ 

au{ A) art (A) 

t e \\ _ a ^(A) _ a r 2 (A) 
[) a u (X) a rl (A) ’ 

ft(\) = Qr3 ( A ) _ a ^(A) 

1 J o rl (A) a e 4 (A) 


where the {a^-(A)} and the {a r j(A)} denote the entries of the transition matrices 
from the left and from the right, respectively. More precisely, 


(4.6) 


where 


a t i( A) =1 e~ lXz + in {z)dz 

Jr+ 

ae 2 (A) = - f e 2lXv u 0 {y)dy - f e lXz + An (z)dz 

Jr Jr 


aes( A) = / e~* Xv vo(y)dy + / e~ lXz + up (z)dz 

Jr Jr 

a u ( A) = 1+ / e iXz <f> up {z)dz, 

J R+ 


(4.7) <f> dn (~) = / u 0 (y)K dll (y,y + z)dy, $ dn ( 2 ) = / u 0 (y)K dll (y, z - y)dy, 


(4.8) 4> up (z) = / vo(y)I\ up (y, y + z)dy, <f> up (z) = / v 0 (y)K up (y, z - y)dy, 
Jr J —00 

and 


on (A) =1+ / e iXz + An {z)dz 

Jr+ 

MA) = / e 2iXy u 0 {y)dy+ / e lXz + dn (z)dz 


(4.9) 


a r 3 (A) = 


= e 


—2i\y 


v o(y)dy — / e~ lXz + up (z)dz 


a r 4 (A) — 1 


* p * 11 

/ e~ iXz ^ np (z)dz 
J R+ 











SCATTERING DATA COMPUTATION FOR THE ZAKHAROV-SHABAT SYSTEM 


7 


(4.10) 

* dn (*) = [ uo(y)M dn (y, y — z)dy, 
JR 


(4.11) 

^ up (2) 


r*+oo 


vo(y)M up (y, z — y)dy, 



u 0 {y)M dn (y,z - y)dy, 


V up (z ) = [ v 0 (y)M up (y, y - z)dy. 

Jr 


While the approximation of T simply requires the computation of at 4 (A) and 
a r i(A), that of p and £ is more complicated. In fact, to approximate p(a) and £(a) 
we first have to compute the scattering coefficients by means of CLl-tEO), then 
the reflection coefficients R( A) and L{ A) by using (14.511 and (14.411 and, finally, p(a) 
and t(a) by resorting to the inverse and direct Fourier transforms as indicated in 
( 12111 ) and (tUTUp . 

The stability of this numerical procedure essentially depends on the decay of 
R( A) and L{ A) for A —> ±00 since the smoother the initial potential the faster their 
decay. If the initial potential has jump discontinuities then R and L decay as A -1 
for A —> 00 while if u 0 £ C°° (M) then R and L decay superpolynomially. 

Hence, this procedure is effective whenever the initial potential is smooth enough, 
that is at least u 0 £ C( R). If this is not the case the Fourier transforms p(a) 
and £(a) could be approximated by solving structured Fredholm integral equations 
stated in the following theorems. The development of an effective algorithm for 
solving these equations is devoted to a subsequent paper. 


Theorem 4.1. The function p{a) is the unique solution of each of the following 
Fredholm integral equation of the second kind : 

n OO -1 

(4.12) p(a) + J o $ up (z) p{z + a)dz = -~v 0 (|) - $“ p (a), 

(4.13) p(a)+y ^ dn (z) p(z + a)dz = ~v 0 - ^ up {a), 

where Q up and are given in (I4.7II - II4.8H and if* 1 and 'b up are defined in (14.101) - 

grlD . 

Proof. Let us first note that from (SSD 


(4.14) a ei (\)R{\) = -a«( A) 

where au and at 3 are dehned in (14.61) . Introducing the Heaviside function H{z) = 1 
for 2 > 0 and H{z) = 0 for z < 0, taking into account that 

i?(A)= R{p{a)} 

and using (14.141) we can write 

(1 + _F { $up (- a )H(-a)})F{p(a)} = -F 0 (|) + $ up (a)| . 

Hence, applying the inverse Fourier transform and the convolution theorem, we 
have 

p(a) + (4> up (-a)F7(-a)) * p(a) = v 0 - $ up (o), 

and then the equation (14.121) is an immediate consequence of the convolution defini¬ 
tion and the Heaviside function. Equation (14.131) can be obtained similarly, noting 
that R( A) satisfies the relation 


Hrl(A)7?(A) — U r 3(A) 
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and that 

a r i(A) = 1 + J r {4 ,dn (-a)iJ(-a)} and a r3 (A) = T j^vo + ^ up (a)j . 

□ 

We note that, from the numerical point of view, it is irrelevant if we solve (14.121) 
rather than (14.1311 . since both are Fredholm integral equations of the second kind, 
equally structured. 

Applying the same technique we obtain the analogous 

Theorem 4.2. The function £{a) is the unique solution of the two structured Fred¬ 
holm integral equations of the second kind: 

(4.15) £(a) + f $ up (z)£(z + a)dz =—^-uo — & dn (a) 

Jo 2 V 2 / 

(4.16) + / * dn (z)£(z + a)dz = ~u 0 (^)-* dn (a), 

where ’F dn is defined in (14.101) and <F“ P and are given in (14.71) - (14.81) . 

We omit the proof, as it is analogous to the previous one, after noting that 
L(X) = T{£(-a)}. 


5. The numerical method 

Let us now assume, for computational simplicity, that the support of the initial 
solution is bounded, that is 

(5.1) uo{x) = 0, for \x\ > L , 

which can be considered acceptable whenever uq(x) —> 0 for |x| —> oo, provided 
that L is taken large enough. This hypothesis, as in part already proved in [Tj, 
allows us to greatly simplify the algorithms for the computation of the auxiliary 
functions and also those for the computation of the Marchenko kernels and the 
Fourier transforms of the reflection coefficients. 

The method we propose provides successively the numerical solution of: 

(1) the four systems (13.51) - (13.81) of Volterra integral equations for the computa¬ 
tion of the four pairs of auxiliary functions; 

(2) the two Volterra integral equations (14.1I) - (I4.2I) for the computation of the 
Marchenko kernels from the left and from the right fand Q r , respectively; 

(3) the transition matrices from the left and from the right, the scattering 
matrix and then the inverse Fourier transforms p of the reflection coefficients 
from the right R and the Fourier transform £ of the reflection coefficient 
from the left L. 

Once the Marchenko kernels f \(a) and f l r (a) and the functions p(a) and £(a) 
have been obtained, the bound states {Aj}” =1 with their multiplicities {mj}™ =1 and 
the norming constants {(F^)y s , (r,.)^} are computed by applying to the monomial- 
exponential sums (I2.5I) - (I2.6I) the matrix-pencil method proposed in Jj and [ET. 

5.1. Auxiliary functions computation. As said before, our numerical method 
for the solution of the Volterra systems m-m is greatly influenced by the 
hypothesis (ELD- It implies a reduction of the auxiliary function supports, which 
allows us to develop algorithms that are simpler and numerically stable. 

As proved in [7], K up and K dn have the supports depicted in Figure [L] Taking 
into account the symmetry properties (13.131) or (13.141) of systems (13.51) and (13.61) , it 
is immediate to check that supp(K up ) = supp(K dn ) and supp(K dn ) = supp(K up ). 
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FIGURE 1. Supports of the auxiliary functions A' up and A' dn (to the 
left) and A' dn and A' up (to the right) 



Figure 2. Additional properties of A' up and A' dn (to the left) and 
K dn and K up (to the right) 


For the numerical solution of system (EH), the following properties, proved in 
[7], are also important: 

1. If x < —A, whatever h, K up (x,y) and K dn (x,y ) are both constant on 
the line y = x + h. For this reason we put K up (x,x + h) = C^, and 
A' dn (x, x + h) = C^ h for each given value h. 

2. If x < ~L and x + y > —2A, K dn (x, y) and I\ up are both constant on each 
line x + y = —2(A — h) for each 0 < h < 2 L. 

These two results are graphically represented in Figure EJ where K dn (x,x + h) = 
C^ h and K up (x, x + h) = C^ h . 

Analogous considerations, based on results reported in [71 allow us to claim that 
the supports of (M up ,M dn ) are those depicted in Figure [5] As for (Ji up , A" dn ) and 
(AT up ,A' dn ) as for the pairs (Af up ,M dn ) we have additional properties very useful 
from the numerical point of view. With obvious meaning of the symbols, they are 
reported in Figure [I] 

A simple inspection of Figures |T] and [3] makes it evident that the area where we 
need to compute I\ up and K dn , as well as K up and AT dn , is given by the orange 
triangle represented in Figure [5] In the remaining areas of the respective supports 
their values are immediately obtained by using those of the orange triangle. The 
orange line shows, in particular, the values of the orange triangle we use to compute 
(.A up ,A' dn ) and (AT up ,A' dn ) in the point of the gray area. Similar considerations 
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FIGURE 3. Supports of the auxiliary functions M up and M du (to the 
left) and M dn and M up (to the right) 



Figure 4. Additional properties of M up and M dn (to the left) and 
M dn and M up (to the right) 



Figure 5. Geometrical visualization of the computational area of 
K up and K du (to the left) and K up and K An (to the right) 

hold true for the computational area of the pairs (M up , M dn ) and ( M up , Af dn ) which 
is depicted in Figure [Gj with the analogous meaning of the symbols. 

Algorithm. Given the initial solution u o and Vo = Uq in the focusing case or 
vo = — Uq in the defocusing case, we have to solve Volterra systems m-m • 
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Figure 6 . Geometrical visualization of the computational area of 
M up and M dn (to the left) and M up and M dn (to the right) 



Figure 7. Sorting visualization of collocation points in the trian¬ 
gle of A' up , A" dn , K up and K dn 



Figure 8. Sorting visualization of collocation points in the trangle 
of M up , M dn , M up and M dn 


Let us start with the numerical solution of system m - As noted before, under 
the hypothesis eu), we can limit ourselves to solve this system in the triangular 
computational area represented in Figure 0 as the values of K up and K dn in the 
remaining parts of their support are then automatically known. 

The algorithm that we propose in this paper is more effective that the one re¬ 
ported in |7 whose aim was simply to check the effectiveness of our approach, 























12 


L. FERMO, C. VAN DER MEE AND S. SEATZU 


highlighting the mathematical problems to overcome to obtain a satisfactory solu¬ 
tion of the problem. Though the collocation strategy is the same used in [7], the 
algorithm used here is more complex and effective. In fact, it is based on the com¬ 
bined use of the trapezoidal rule, the composite Simpson quadrature formula and 
the 3/8 Simpson quadrature rule JT2J Section 3.1], instead of only the composite 
trapezoidal quadrature formula used there. 

The first step is to fix a proper mesh in the computational area which can be 
done by fixing ngN, taking h = ^ and introducing the following mesh points: 

£>fc = {{xi,x i+ 2k), Xi=ih, i = n - k,n - k - 1,..., —n + 1, -n} 


where the index k = 0 , ..., 2 n identifies the line y = x + 2 kh on which we want 
to compute the unknown functions, whereas i labels the abscissa of the i- th mesh 
point on the line. 

For the sake of simplicity, let us hereafter write u and v in place of uq and i'o, 
respectively. The computational strategy requires us to compute first K up and 
A dn in the nodal points of the bisector (xi,Xi). Consequently, recalling (E2D and 
denoting by A'//, , K dn a the approximation of K up (x, y), K dn (x , y) in the nodal points 
of T> o, we can write 

i r°o i rx n+ i 

K Ti=-^J u(z)v(z) dz = -- j u(z)v(z) dz 

Kt% = -jVu i = n,n — 1, ..., —n + 1, —n. 

To approximate the above integral, it is convenient to use different quadrature 
formulae, according to the node x t . More precisely for: 

• i = n, being involved only two nodal points, we use the trapezoidal rule 

A'/i = ~{u n v n +u n+ iv n+1 } = ~u n v n 
as, for (13.9[) . u n+ iv n+ i = 0; 

• i = n — t, i = 1, 3, 5, ... , 2n — 1, we apply the composite Simpson rule. 
Recalling that u n +i = u n +i = 0, we then obtain 


K up = — 

*.* 3 


~r J- _1 

2 2 1 

UjVj + 4 Uj+2j-lVj+2j-l + 2 Ui+2jVi+2j 

.7=1 .7=1 


• i = n — £, £ = 2, 4, 6 , ..., 2n, noting that 




u(z)v(z) dz 



[■x n +i 

A i+3 


u(z)v(z)dz, 


and that the first integral involves four nodes, while the second involves an 
odd number of nodes, we can apply the 3/8 Simpson rule [T2] p.128] for 
computing the first integral and the composite Simpson quadrature formula 
for the second one. Hence, recalling again that u n + 1 = v n +i = 0, we have 

3 

K i P i = o h l UiVi + 3u.j + iU i+ i + 3u l+2 V i+ 2 + Ui +3 Vi +3 } 

’ O 

L L -1 1 

2 2 x 

Ui+ 3 Vi+ 3 + 4^^ Ui+2+2j + 2 Ui+ 3 +2jVi+ 3 +2j ■ 

.7=1 i=i 



Once A up and K dn on the nodal points of the bisector y = x are known, to 
evaluate them on the nodal points of the parallel lines to the bisector, we collocate 
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system (E3D on the nodes of the mesh (xj, Xi+ 2 k), taking successively k = 1 , ..., 2n 
and, fixing k , assuming * = n — k, ..., — n + 1, — n. Hence, we can write 

/•OO 

K pp i+ 2 k + / U (X) K du (z, z + 2kh) dz = 0, 


K dn 

I *’i,i+2k 


v(z) K up (z, 2(i + k)h — z) dz = \vi + k- 


These formulae, taking into account the support of the functions involved (Figure 
01 , reduce to 

f'Xn — fc + 1 

u{z) k dn (z, z + 2 kh) dz = 0 , 


K up 


K dn 


;(z) K up (z,2(i + k)h - z) dz = \v i+ k- 


To compute the first integral 

rXn-k+1 

iL = 


u(z) K dn (z, z + 2kh) dz 


we use different quadrature formulae, according to the node X{. More precisely, 
fixing fc, for: 

• i = n — k, being involved only two nodal points, we use the trapezoidal rule 
and then take 

1k,i t ~2^ri—k^n — k,n+k^ 

as the nodal point (x n -k+i, x n +k+i) is outside of the support of Jv dn (x,y); 

• i — n — k — £, with £ odd and £ < 2n — fc, applying the composite Simpson’s 
rule, we obtain 

i+i 


ri - 

1 k,i — 


u i^i^i+2k + u i+2j-lKi+2j-l,i+2j-l+2k 
.7=1 


>dn 


as K dn _ k+l n+k+ i 


Ui+ 2i_i A, 
= 0 . 


dn 


+2 u i+2j-l-t\ i+2 j.i+2j+2k 

j= 1 


■i = n — k — £, with £ even and £ < 2 n — k 1 noting that 


H _ 

1 k,i ~ 


f 

J X r 


u(z) K dn (z , 2 + 2kh) dz , 


we apply the 3/8 Simpson’s rule for the first integral and the composite 
Simpson’s quadrature formula for the second one. Hence, we have 
3 - - 

= O ^ \ u i^i,i+2k + ^ u i+l-^i+l,i+l+2k + 3 w i+2-?Q+2,i+2+2fc + M i+3^j+3,i+3+2fc] 


Ui+3^+3 ,i+3+2fc + 4 El M »+2+2 i i^'u+2+2?.i+2+27+2fc 
1 = ! 

e 


4-1 


+2 E ^+3+2j-^f+3+2j, 


i+3+2j+2k 










?s- to 
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as the nodal point (x n -k+i, x n +k+i) is outside the support of K dn (x,y). 
The computation of the second integral 


fXi+k 


T- — 

1 k,i — 


v(z) K up (z, 2(i + k)h 


z) dz, 


is also based on the use of quadrature formulae, essentially dependent on the line 
y = x + 2kh. More precisely, for: 




k = 1 , as only two nodal points are involved, we apply the trapezoidal rule, 
obtaining 


^ k,i ~ 2 + w *+l-Kj+l,*+l}) 


• k = 2,4, 6 ,2 n, we use the composite Simpson quadrature formula. Pro¬ 
ceeding in this way we obtain for i = n — k ,..., — n 


L 
2 

V ^i^i+2k + l,i—2j+l+2fc 

j=l 

1 

+2 E v i+2jKi+2j,i-2j+2k + V i+k^i+k,i+k 

j = 1 

= 3 V *^>i+2fc + w k,ii 

where Wky is the sum of the K up values in the nodal points belonging to 
the bisector and the previous parallels. In fact, the K up values of the first 
sum belong to the lines y = x + [2k — 2(2j — 1 )\h, those of the second one 
belong to the lines y = x + [2k — 4 j]h and the last term to y = x. 

• k = 3,5,7,..., 2n — 1, we write 



T 2 = 
k.i 



v(z) K up (z, 2(i + k)h — z) dz 


and then we use the 3/8 Simpson rule for the first integral and again the 
composite Simpson quadrature formula for the second one: 


= -h 

8 

h 

+ 3 


,i+2fc-l + '^ Vi + 2 ^i+2,i+2k-2 + V i+^i+3, i+2fc-3 


S+3,i+2fc-3 + ^y, V i+2+2j K d 2+ 2i.i-2-2 1+ 2k 
3 =1 


+2 ^ u i+3+2jK i+3+ 2j t i-3~2j+2k + Vi + k ^i+k,i+k 

3 =1 

3 

= gf lv i^iy+2k + w k>*> 


where Wk,i is known, being a linear combination of -ftT up values already 
computed. 

Once the integrals have been approximated as described above, we obtain 
the 2 n following structured systems of order 2 ( 2 n + 1 — k), k = 1 , ..., 2 n 


K p +_ u M kf = o 

U fc , 2 k/ P + k/ n = v fc - w fc 


(5.2) 
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that allow us to compute the functions K up and K dn in the 2n +1 — k nodal 
points of T>k as 


iup _ /T^up 

— v fY n-/s,n+A;> 

f/dn _ / dn 

y-^n—k,n+ki 


K up 

v n— k— l,n+k —1 ’ 

K dn 

x v n— k—l,n+k— 1 ’ 


K up K up ) T 

• • ’ ^ -n+l,-ri+2fc+l> -n,-n+2k) 

TS U P ff dn 

■ ' > - fi -n+l-ra+2fc+l> xy -n,-n+2kJ • 


Notice that Ufc,i, Ufc ,2 are the following structured matrices: 
Ufc ; 2 = Ck hdiag(v n -k, Vn—k—ii ..., v — n ) 


with ci = 1/2, C 2 = C 4 = ... = C 2 n = 1/3 and 03 = 05 = ... = C 2 n -i = 3/8 
and 


/ 1/2 




4/3 

1/3 



9/8 

9/8 

3/8 


4/3 

2/3 

4/3 

1/3 

4/3 

17/24 

9/8 

9/8 


k, 1 — h 


4/3 2/3 . 17/24 9/8 9/8 3/8 

\4/3 2/3 4/3 . 4/3 2/3 4/3 1/3/ 

diag(u n -k, Un—k — 1: • * • ) n+1: ^—n)' 

The most obvious computational strategy is to reduce to a sequence 
of n — k systems of order two. However, our numerical experiments indicate 
that the numerical stability increases by using a suitable iterative method. 
It requires solving iteratively the system 

(5.3) (I-U M U M )iq p = U M w fe , 

and then computing 

(5-4) K/ n = V fc —U M K" p . 

The matrix of system (15.31) . for h small enough, is diagonally dominant as 
each nonzero element of contains a factor h 2 , so that the Gauss- 

Seidel method is a suitable choice of iteration method, assuming as an 
initial vector the values of k/ p in the previous parallel, that is taking in the 
(k + l)th parallel to the bisector 

(5.5) (k^!) (0 ) =k pp fc = 0,l, ...,2n-l. 

As I — is lower triangular, it is of course possible to solve it by a 

descending technique. 

Remark 3. Once we have solved system m we can immediately deduce the so¬ 
lution of system (EH) taking into account Remark [1] In any case, we note that, as 
the computational area of system EH is the same as that of EH , the algorithm 
to solve (13.61) is analogous to that adopted for system (13.51) . 

The same comparative considerations hold true for the computation of (M up , 
M dn ) and (M up , M dn ) in the nodal points of their computational area. Moreover, 
although the computational area for (M up , M dn ) is not the same as that for ( K up , 
A' dn ), the technique for their computation is essentially the same. 

Noting that (Figures [5j H the two computational areas are symmetric with 
respect to each other, we first have to compute (M up , M dn ) in the bisector and 
then on the parallel lines y = x — 2kh, k = 1,2,..., 2 n. Furthermore, to compute 
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M dn in the bisector we can adopt the same algorithm for K up as relations (13.911 and 
(13.1211 indicate. A comparison between the systems (13.51) and (13.71) also suggests 
to approximate the first integral in (ITTb by a simple adaptation of the method 
developed for the second one in (13.51) . as well as the second integral of (13.71) by 
adapting the method for the first integral of (13.51) . 

5.2. Marchenko kernel computation. To compute and O r , that is to solve 
the integral equations 4 HD and (El) . we first note that ED implies the bounded¬ 
ness of their supports. In fact, as proved in [TJ Lemma 5. i], ED implies that 

suppiyit.) = [0, 2 L\, and supp[VL r ) = [—2L, 0]. 

For the approximation of we collocate ED in the nodal points 

{ (%n—2i 3 %n) ; %n—2i (jl 2i)/l, i — 0,1,..., 77.}, 


by obtaining 

(5.6) (^'2(n—i )) 3“ / (x n —2ii Z^Qg^Z — K (^(n— 2i) 5 *£n)- 

J X n -2i 

Now, to compute the above integral we use different quadrature formula by adopting 
a steplenght 5 = 2h that is twice the one considered in the numerical solution of 
system ED to avoid the interpolation among the values of the auxiliary functions 
computed before. More precisely, for 

• i = 0 we immediately obtain that 

= ~K nn = — — V n 


in virtue of (13.91) : 

• i = lwe use the trapezoidal rule by getting 


1 + ^K-2,n-2 ) fy,2(„-l) = -K-2,n ~ K-2^2n, 


• i = 2,4,6,... we use the Simpson quadrature formula 


1 4 - — h' dn 
1 ' ^ Iv n—2i,n—2i 


Q{,2(n—i) — K n -2i,n g 1^ 0 ^ ^n-2i.n—2(i- ?) ^i,2(n-2(i—j)) 

3 =1 


+ 2 X] K n-2i,n-2(i-j-l)^e,2{n-(i-j-l)) + I^n-2i,n^(,2n 

3=1 

i = 3,5, 7,... as we can write 

/ K dn (x n - 2 i,z)fle(z + x n )dz 


r X "-2(i-3) 


■L 


K dn (x n - 2 i, z)ttt(z + x n )dz 


x n — 2(i — 3) . 


we approximate the first integral by using the 3/8 Simpson rule and the last 
integral by adopting the composite Simpson quadrature formula. Hence we 
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get 


3<5 


3(5 


1 d“ ^ J ^£,2(n—i) 2i,n g i^^n—2i,n—2{i— l)H^,2(7i—(i—1)) 

+ ^ K n-2i,n-2(i-2)^e,2(n-(i-2)) + -^n-2i,n-2(i-3)^,2(n-(i-3))) 


- - K 


n—2i,n—2(i 


-3)^,2(n-(i-3)) + ^n-2i,n-2(2i-3-j)^,2(n-(2i-3- 


1=1 


!--l 


+2 -^n-2i,ri-2(2i-4-j)^,2(n-(2i-4-j)) + ^n-2i,n^i,2n 

1=1 

An analogous procedure can be applied to approximate S2 r in [— 2L, 0]. More 
precisely, we collocate (14.21) in the nodal points 

{(x 2 i-„,a;-n), x 2 i~ n = (2i - n)h, i = 0, 1, n}, 

by obtaining 

r x 2 i-n 

( 5 . 7 ) tt r (x 2 (i-n))+ M up (x 2 i-n,z)fle(z + X- n )dz = -M up (x(2i-n),X-n)- 

J X-n 

Hence, by adopting the technique illustrated above, 

• for i = 0 we immediately obtain 

Hr,—2n = = ——U-n 

in virtue of fl3H) ; 

• for i = 1 we obtain 


1+ 2 ^2-n,2-nJ H r ,2( 1 _n) — ^2-n,-n ^2-n,-rS^r,-2n\ 

• for i = 2,4, 6,... we obtain 

+ ^^2i-n,2i-n^ ^r,2(i-n) = ~^2i-n,-n ~~ 3 ^4 ^ ^2i-n.2(i— i' 1 -n^r, 2 ((»-j)-r 

i±i-i N 

+2 ^ ^2i-n,2{i-j-l)-n^r,2({i-j-l)-n) + M2i-n,-n^r-2n 


1 = 1 


• for * = 3,5,7,... as we can write 



M up (x 2 i-n, Z)tt r (z + X- n )dz 


{‘ x 2(i — 3) — n 



M up (x 2 i -n,z)n r (z + X- n )dz, 


we approximate the first integral by using the composite Simpson rule and 
the second one by adopting the 3/8 Simpson’s quadrature formula. Hence 


i)) 
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we get 


1 + ^M up 


2i-n,2i-n J n) — lv± 2i—n,—n 

UP ^r.2((i-2)-n) + Mf? _ o/.- « 


= -M“P 


35 / 

"8" V 


3-^2i > -n,2(i_l)_n^D2((i-l)-n) 


+3M, 

_ 5 

~~ 3 


.) 


2i-n,2(i-2)-n Lr ,2(.( i -2)-n) T lvl 2i-n,2(i-3)-n V,2((i-3)-? 

t+i 

^2i-n,2(i-3)-n^r,2((i-3)-n) + 4 ^ P _ 3 2 f 2 i-3—il-n ^r,2((2i-3-.7')-n) 

J=1 


+2 E M-n— 2i,2(2i-4-j)-rS t ‘r,2((2i-4-j)-n) + ^ 2 i-n,-n^,-2n ) • 

j=l 


5.3. Computation of the scattering matrix and inverse Fourier transforms 
of reflection coefficients. In this section we illustrate our method to approximate 
the scattering matrix and then to compute the transmission coefficients T defined in 
(14.31) . the reflection coefficients R and L introduced in (I4.5I) - (I4.4I) and their Fourier 
transforms p and t given in (I2.9I) - (12.10I) . under the assumption that uq € C(R). 


Approximation of the transmission coefficient T. It is based on the two 

equivalent definitions of the transmission coefficient 


(5.8) 


n a) = 


i 


om( A) ’ 


T{ A) = 


1 


a r i(A) 


that is on the computation of the coefficients of the transition matrices 

(5.9) a u (\) = l+ [ e iXz <fr up (z)dz = 1 + 2itJ 7 ~ 1 {$ up (A)Kf(A)} , 

(5.10) a r i(A) = 1 + f e iXz ty dn (z)dz = 1 + 2irJ 7 ~ 1 {^ dn (X)H(X)} , 

J R+ 

where H denotes the Heaviside function and J 7-1 {g} stands for the inverse Fourier 
transform of g. 

Let us only illustrated the algorithm for the computation of the coefficient an 4 
as the computation of a r 1 is analogous. 

At first we note that, taking into account (15.11) and the support of K up , the 
kernel <f> up of (14.12|) can be written as follows: 


(5.11) 


<S> up (z) 



vo{y)K up {y,y + z)dz, 


0, 


for 0 < z < 4L 
for z > 4L. 


Then, writing, for simplicity, <J>“ P = *l> up (^) = <!> up (2 hj), j = 0,1, 2,..., 2n— 1 we 
have successively to compute <E>q P , $ pp , ..., $2n-i by obtaining 

/ (■ n-j)h 

v 0 (y)K pp (y , y + 2 hj)dy, j = 0,1,..., 2 n - 1. 

-nh 

We remark that its computation requires only the values of A' up (y, y + 2 hj) which 
we have already computed since they are the values of K up on the jth parallel to 
the bisector y = x. For this reason $“ p can be computed by simply adopting the 
computational strategy that we developed for computing K up . At this point the 
approximation of T(A), easily follows by using (15.81) . 
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Approximation of the reflection coefficients R and L. In the matter of the 
computation of the reflection coefficients, taking into account flUD and (TOD , we 
can write 


(5.12) R(A) = -T(A)o, s (A), L(A) = T(A) o n (A) 

where T(A) = ^y 

o«(A) = [ e~ iXy 

Jr 

and 

002 (A) = ~ f R elXV (f) + $dn (y)^ d V = j^Mo (|) +^ dn (2/)| • 

Other equivalent expressions can be deducted by using the definitions of R , L and 
T in terms of the coefficients of the transmission matrix from the right. 

To approximate <^ 3 , taking into account (EH and the support of K up , first we 
note that 


Vo (f) + ^ UP(2/) 


dy = T 


-vo 


(|) + $ up (y) 


(5.13) 


<$> up (z) = 


Vo{y)K up (y, z - y)dy, for \z\ < 2 L 
for \z\ > 2 L. 


Moreover, adopting the notation used before and noting that = <b up (— n) = 
i ,up (— 2 nh) = 0 , we can write 

= f v 0 (y)K up (y,2hi - y)dy, i = -n + 1,..., 0, ..., n. 

J —nh 


Hence p , as well as $“ p , can be computed by simply adapting the computational 
strategy developed for K up . The approximation of R and L immediately follow by 
using (15.1211 . 


5.4. Computation of the bound states and the norming constants. For the 

sake of completeness, we now give a brief description of the matrix-pencil method 
that we have recently developed for the identification of the bound states and 
the norming constants [ 610 . Setting Zj = e lXj , the spectral function sum Se(a) 
introduced in (12.51) can be represented as the monomial-power sum 

n mj — l 

St(a) c Js a s z°, 0° = 1. 

j =1 s=0 

Letting M = m\ + ... + m n , the method allows one to compute the parameters 
{n,rrij,Zj} and the coefficients {c JS }, given Se(a) in 2 N integer values (N > M) 

a = ao, cko + 1, • ■ • 1 cxo + 2 N — 1 , with ao G N + = {0,1, 2,...}, 

under the assumption that a reasonable overestimate of M is known. 

The basic idea of the method is the interpretation of Se(a) as the general solution 
of a homogeneous linear difference equation of order M 

M 

^ ^ PkRk+ap = 0 
k—0 

whose characteristic polynomial (Prony’s polynomial) 

n M 

p ( ,) = n(-^r = E^ P" si 

j=1 k —0 
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is uniquely characterized by the Zj values we are looking for. The identification 
of the zeros {zj} allows one to compute the coefficients Cj s by solving in the least 
squares sense a linear system. 

For the computation of {zj} and then of the bound states A j, the given data are 
arranged in the two Hankel matrices of order N 

(S° t )ij = S t (i+j-2), (S 1 e ) ij = S e (i + j-l), i,j = 1,2,..., N. 


To these matrices we then associate the M x M matrix-pencil 
S mm(z) = (Sjvm — z &nm) 


where the asterisk denotes the conjugate transpose. As proved in [5], the zeros 
Zj of the Prony polynomial, with their multiplicities, are exactly the generalized 
eigenvalues of the matrix-pencil Smm(~)- The simultaneous factorization of the 
matrices S° NM and S l NM by the Generalized Singular Value Decomposition allows 
us to compute the zeros Zj and then the bound states A j, as A j = —ilogZj. 

Analogous results can be obtained by a proper factorization of the augmented 
Hankel matrix = [S° x , S|], where S° 1 is the first column of S^ M and S° is 
obtained by S* by simply deleting its last column. As shown in IB], the QR fac¬ 
torization of S e is as effective as its SVD factorization considered in jS], though its 
computational complexity is generally smaller. 

The vector of coefficients 


C — [Cl 0j ni — 1, *■•; C-L 0; CL fiM-l] 

is then computed by solving (in the least square sense) the overdetermined linear 
system 

X2-0 „ __ c0 

where S^° = [S^(0), Se(l), ..., Sg(N — 1)] T and K^ M is the Casorati matrix asso¬ 
ciated to the monomial powers {k s Zj} for k = 1, ...,N — 1. 

If rrij = 1, the Casorati matrix reduces to the Vandermonde matrix ( V)ij = 
Zj +1 of order N xn associated to the zeros z i, ..., z n . The solution of the Casorati 
system allows us to immediately compute the norming constants as (r^s = s!c Js . 

The coefficients {(T r )j S } are then obtained by solving, in the least square sense, 
a linear system whose vector of known data is given by Q r (a) evaluated in a set of 
N points, with a sufficiently large N > M. 


6. Examples 

Let us now present two examples. The first one is a reflectionless case while the 
second one has reflection coefficients different from zero. Each of them will be used 
in the next section to give a numerical evidence of the effectiveness of our method. 

Example 1 (One soliton potential) 

Considering the initial potential for the NLS in the focusing case we take 

(6.1) Uo(x) = 2ii r ie 1 ^ x+ ^ sech(xo — 2 tix) 

where £, 4>, Xo £ R and 0 ^ r] £ R. As proved in [9], the corresponding initial value 
problem (ED can be solved exactly, as already considered in several papers and in 
particular in [5] and [1|. Let us note that 277 > 0 represents the amplitude of the 
initial potential and /io = Xo/2r] is the initial peak position. 

In this example the norming constants from the left and from the right are 0]: 

(6.2) T^ = 2irie x °-^ and T r = -2ir]e- Xo+i ^. 
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Moreover, setting a = 77 + i£, as it is immediate to check, the exact solution of 
the Volterra system EH for y > x is 

_p* e ~ a *( x +v) 


^' up ( Z, yp 

yK dn (x,y), 


l 4 . £ 2 (x 0 - 2 r/x) 


(x+y)—2ax ] J 

2r/ 


while the exact solution of system EH can be obtained by resorting to relation 
(13.1311 . Furthermore, the closed form solution of the Volterra system (13.811 is 

'|iv 2 


<M up (x,yy 


1 


1 4 - e -2(x 0 -2rix) 


0 a* (x+y)+2ax 

2 V " 

_r* e a *( x +v ) 


^M dn (x,y); 

while the solution of system (13.71) can be deducted by using relation (13.131) . As it 
represents a reflectionless case, 

p(a) = 1(a) = 0, a € R, 
so that the exact initial Marchenko kernels are [4] 

tte(x) = r t e~ 
f l r (x) = r r e Q 

Finally, the scattering matrix is 

( A + ia* 

S (\\={ T W 


0 —ax 
" 1 

„ax 


\R(X) T(X)J 


A — ia* 

0 


0 

A + ia* 


A e C+. 


A — ia*, 

Example 2 (Gaussian potential) 

As a second example of the initial potential for the NLS, we take 

(6.3) u 0 (x)=q 0 e^ x e~^, 

where go > 0, a > 0 and 

As in mm we investigate the defocusing case in which the scattering coeffi¬ 
cients T(A), R( A) and L( A) are all continuous functions and there are no bound 
states. Hence, in this case the following relations hold true 


(6.4) 


£lg(a) = p(a) f l r (a) = i(a). 


qo^/ncr < 


Moreover, we also consider the focusing case. In such a case, whenever 

7r 
2 ’ 

there are no discrete eigenvalues. On the contrary we have n discrete eigenvalues, 
all of them simple and having real part — if m 


(6.5) 


1 


n — — 7T < goV 7rcr < I n + n I 7T- 


1 


As a result the spectral sums from the left and from the right (12.51) and (12.61) reduce 


to 


(6.6) 

n 

St(a) = 

j=i 

(6.7) 

s r (a) = ^(r r ) je iA 7*« 


i=1 


a < 0. 
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We also remark that the reflection coefficients R( A) and L( A) and the transmission 
coefficient T( A) are discontinuous at A = — ^ if [TU] 


9ov 7rcr = 



7T 


for some positive integer n. 


7. Numerical results and conclusions 


Test 1 (One soliton potential) 


Let us consider as in [?] the initial potential (16.11) with £ = 1/10, xq = <j> = 0 
and r) = 2. In order to compute the non-zero scattering parameters that in this 
case are the norming constants, the bound states and the transmission coefficient, 
at first we solve the Volterra’s system (13.61) and (13.81) with L = 8 and n = 3000 by 
obtaining the following relative errors 


||A' up - AT up || 

pwi 

\\M up - M up || 

ll^ up ll 


= 1.80e — 06, 


= 1.07e- 07, 


|| K dn - /\ dn || 
|| M dn - M dn \\ 


1.04e - 07, 


1.80e- 06, 


where here and in the sequel the ~ sign denotes the approximation of the exact 
function previously given and || • || denotes the maximum norm of the involved func¬ 
tion in their computational areas. Identical relative errors are of course obtained 
for the remainding auxiliary functions, as a result of the symmetry properties (13.131) 
and (13.141) . 

Once these auxiliary functions are computed we numerically solve equations 
getting for the Marchenko kernels from the right and from the left with the following 
relative errors: 


max | VlAx) — Tig (rc) | 

rce[0,2 L] 


max | fig (a:) | 

®G[0,2 L]' 


max |fl r (a:) — fl r (a:)| 

x£[-2L,0] 


max |fL(:r)| 

rcG[-2L,0] 


~ 3.24e - 07, 


where the symbol — means that the left term coincide with the right term up to 
the third decimal digit. 

At this point, by using such kernels, we apply our matrix pencil method 0 by 
finding a single bound state term, a norming constant from the left and a norming 
constant from the right with the following relative errors: 


|A~A| 

|A| 


= 4.lie — 09, 


| ft - T e | 

|r*| 


|r~r-r r | 

IP| 


~ 3.24e- 07. 


In the matter of the relative errors of the scattering matrix, we obtain 


max || S(A) 

Ae[—2L,2L] 


S(A)|| 


max 

Ae[-2L,2L] 


l|S(A)|| 


= 4.60e — 07. 


Moreover to ascertain the effectiveness of our numerical method we checked the 
numerical validity of the algebraic property IE2D- The results are at all satisfactory 
as Figure [9] shows where the behavior of the error function 


E s ( A) 


; l -(S t (A)JS(A) + S(A)JS t (A)) - J 


is reported for A € [—2L, 2 L\ in semilog scale. 
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Figure 9. E s { A) in semi logarithmic scale 


Concerning the trasmission coefficient, we can compute it by approximating at 
first the integral $ up defined in (14.SI) , and then using (15.81) . In Table [T] we give the 
following relative errors we obtain for such a coefficient over segment of width 4L 
of three different lines 


E r (T) 


max |f (A) - T(A)| 

AE[a,6] 


max |T(A)| 

A(E[a,6] 


Table 1. E r (T) in the one soliton case 


[a,b] 

E r (T) 

[-2L, 2L] 

3.13e — 07 

[-2L + i, 2L + i] 

2.21e — 07 

[-2L + 5i,2L + 5i] 

4.47e - 07 


Test 2 (Gaussian potential) 

Let us consider first the initial potential (16.31) in the defocusing case with = 
1.9, n = 1, a = 2 as in OH El- To this end, we compute the solution of sys¬ 
tems (1331)-(IXH1) considering as in the soliton case L = 8 and n = 3000, then we 
solve equations sm-dia), compute the scattering matrix and thus the Fourier 
transforms of the reflection coefficients. Our numerical method recognizes that, as 
theoretically expected, there are no bound states and relations (16.41) are numerically 
satisfied since we have the following errors: 


max IfiAa;) — p(x)\ = 1.08e — 10, max — £(x)\ = 1.44e — 09. 

rce[0,2L] xe[-2L,0] 


As in the one soliton case, we checked if our numerical results satisfy the algebraic 
property CO} for the scattering matrix, by considering in semi logarithmic scale 
the error function 


Egd{ A) 


|(St(A)S(A) + S(A)S+(A)) - I 


for A € [— 2L,2L\. As shown in Figure [TO] its numerical validity is satisfactory as in 
the soliton case. 
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Figure 10. Eqd (to the left) and Eqf (to the right) in semi 
logarithmic scale 


Now let us investigate on the focusing case considering the initial potential (16.31) 
with qo = 2.5, \i = 1, cr = 2. As a result, inequality (16.51) implies that we have 
two simple bound states {Ai,A 2} whose real part is —1/2. At first we compute 
the auxiliary functions by solving systems (I3.5I) - (13.8I) with L = 8 and n = 3000, 
then we solve equations ED-(S3, compute the scattering matrix and the Fourier 
transforms of the reflection coefficients. At this point, we apply the matrix pencil 
method described in Section [5.41 assuming that we have not more than five bound 
states. Our method recognize that, as theoretically expected, we have two simple 
bound states having real part equal to — n/2. In fact we get 

Xi = —0.50 + 1.97i A 2 =-0.50 + 0.79i 


with the corresponding norming constants 

r eA = 9.28 - 1.50 10" 8 i r £j2 = 3.74 - 1.7610~ n i 
r r ,i = 9.28 + 1.50 10 _8 i r rj2 = 3.74 + 1.7610“ n i. 

Finally, in Figure [lU] we represent in semi logarithmic scale the error function 


Egf{ A) = 


i(S t (A)JS(A) + S(A)JS t (A)) - J 


for A £ [—2L, 2 L\ that we have computed to check the validity of the algebraic 
property ([23- 


Conclusions 

The numerical results show that our numerical method is effective in both the fo¬ 
cusing and defocusing cases, provided the initial potential decays to zero at infinity 
and is at least continuous. This positive result is due to the possibility to know each 
pair of functions on the whole plane, by solving the relative Volterra system on a 
bounded computational triangle. The accuracy of the identification of the spectral 
parameters strongly depends on this result, since all the subsequent computations 
require the knowledge of the auxiliary functions on their computational triangles. 

We believe that the method can be extended, with the same accuracy of the 
results, in the presence of jump discontinuities of the initial potential. To this 
end, a numerically stable method for the solution of Fredholm integral equations 
(I4.12I) - (I4.13I) and (14.15I) - (I4.16I) should be developed. The development of such a 
method should also be accompanied by an extensive numerical experimentation 
which requires the exact knowledge of scattering data in at least one case in which 
the initial potential has jump discontinuities. Considering that such research takes 
a rather long time, the development of such a method is postponed to a next paper. 
























SCATTERING DATA COMPUTATION FOR THE ZAKHAROV-SHABAT SYSTEM 


25 


8. Appendix 


Supports of the auxiliary functions 

In this section we determine the supports of the auxiliary functions K(x,y) and 
M[x,y) if the potentials uq(x) and vq{x ) have their supports in [—L,L\. It suffices 
to prove parts (2) of Lemmas 5.1 and 5.2 in [7], because the proofs of the other 
three parts of these two lemmas are immediate and proceed as in the discrete case. 
Put 

poo poo 

v(K up \x) = / \K up (x, y)\ dy 7 is(K du ; x)= \K dn (x,y)\dy; 

J X J X 


Q(x) = max(|u 0 (a:)|, |w 0 (a;)|), P(x) = u(AT up ; x) + v(K dn ; x), 

where Q and P are bounded m Then for x < L and x + y > 2 L the integral 
equations (3.1) have zero right-hand sides, because vo(^(x + y)) = 0 for x + y > 
2 L. Integrating the absolute values of K up (x,y) and K dn (x,y) with respect to 
y € (x, +oo), we obtain 


u(/\ up ; x) < / \u 0 (z)\is(K dn -,z)dz, 


v(K dn ] x) < / \v 0 (z)\p(K up -z)dz, 


so that 

0 < P(x) < f Q(z)P(z ) dz. 

J X 

Hence iterating two times the last inequality we have 
pL pL pL pL 

P{x) < / Q(z)P(z)dz< / Q(z) / Q(t) / Q(w)P(w) dw dtdz 

J X j X j Z j t 


< 


Q{w)P{w) dw 


Q(w)P(w) dw 


Q(w)P(w) dw 


Q{z) J Q(t)dtdzj 

J Q(t ) dt j dz 


1 d_ 

2 dz 


Q{t) dt 


z—L 


Q(w)P(w) dw 


Q(t) dt 


Thus iterating n — 1 times we get 


0 < P(x) < 


Q{w) dw 


Q(z)P(z) dz. 


Taking the limit as n —> +oo, we get P(x) = 0 and hence K up (x, y) = K dn (x, y) = 0 
for almost every y > x, as claimed. The proof of part (2) of Lemma 5.2 is analogous. 
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